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Abstract 

The homogeneous isotropic Boltzmann equation (HIBE) is a fundamental dynamic 
model for many applications in thermodynamics, econophysics and sociodynam- 
ics. Despite recent hardware improvements, the solution of the Boltzmann equation 
remains extremely challenging from the computational point of view, in particu- 
lar by deterministic methods (free of stochastic noise). This work aims to improve 
a deterministic direct method recently proposed [V.V. Aristov, Kluwer Academic 
Publishers, 2001] for solving the HIBE with a generic collisional kernel and, in 
particular, for taking care of the late dynamics of the relaxation towards the equi- 
librium. Essentially (a) the original problem is reformulated in terms of particle 
kinetic energy (exact particle number and energy conservation during microscopic 
collisions) and (b) the computation of the relaxation rates is improved by the DVM- 
like correction, where DVM stands for Discrete Velocity Model (ensuring that the 
macroscopic conservation laws are exactly satisfied). Both these corrections make 
possible to derive very accurate reference solutions for this test case. Moreover this 
work aims to distribute an open-source program (called HOMISBDLTZ), which can 
be redistributed and/or modified for dealing with different applications, under the 
terms of the GNU General Public License. The program has been purposely de- 
signed in order to be minimal, not only with regards to the reduced number of lines 
(less than 1,000), but also with regards to the coding style (as simple as possible). 
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The solution is based on the method proposed by Aristov [1], but with two sub- 
stantial improvements: (a) the original problem is reformulated in terms of particle 
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vation during microscopic collisions) and (b) a DVM-like correction (where DVM 
stands for Discrete Velocity Model) is adopted for improving the relaxation rates 
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is particularly important for describing the late dynamics in the relaxation towards 
the equilibrium). Both these corrections make possible to derive very accurate ref- 
erence solutions for this test case. 
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The nonlinear Boltzmann equation is extremely challenging from the computational 
point of view, in particular for deterministic methods, despite the increased com- 
putational power of recent hardware. In this work, only the homogeneous isotropic 
case is considered, for making possible the development of a minimal program (by 
a simple scripting language) and allowing the user to check the advantages of the 
proposed improvements beyond the Aristov's method [1]. The initial conditions are 
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easily modified. 
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From minutes to hours (depending on the adopted discretization of the kinetic en- 
ergy space). For example, on a 64 bit workstation with Intel@ Core^^ i7-820Q 
Quad Core CPU at 1.73 GHz and 8 MBytes of RAM, the provided test run (with 
the corresponding binary data file storing the pre-computed relaxation rates) re- 
quires 154 seconds. 
References: 

[1] V.V. Aristov, Direct Methods for Solving the Boltzmann Equation and Study of 
Nonequilibrium Flows, Kluwer Academic Publishers, 2001. 



1 Introduction 

In a dilute gas, tlie Boltzmann transport equation [T||2] describes the time evo- 
lution of the single-particle distribution function, which provides a statistical 
description about the positions and velocities of the gas molecules. From the 
theoretical point of view, it is one of the most important equations of non- 
equilibrium statistical mechanics and one the most powerful paradigms for ex- 
plaining transport phenomena in fluids. Moreover, from the engineering point 
of view, since early fifties, it received a lot of attention due to aerodynamic 
requirements for high altitude vehicles and vacuum technology requirements 
[2] . Nowadays, the set of applications has been widen by including dilute gas 
flows in micro-electro- mechanical systems (MEMs) [3]. These devices are in- 
creasingly applied to a great variety of industrial and medical problems. In 
these problems, given the small dimensions of the devices, it is necessary to use 
the kinetic theory, instead of the usual fluid dynamics, based on the Navier- 
Stokes equations, to describe the motion of dilute gases in the small gaps of 
these devices. 

Because of the intrinsic complexity of this equation (the single-particle distri- 
bution function is defined in the phase space and the time evolution is ruled 
by a five-fold coUisional integral), solving the nonlinear Boltzmann equation 
is extremely complex. Hence, from the very beginning, there was an attempt 
to formulate simpler models, which preserve the main features of the dynamic 
approach to the thermodynamic equilibrium. As pointed out in the Cercig- 
nani's biographical work [1], Boltzmann himself started in his fundamental 
paper [1] by considering first the case when the distribution function does not 
depend on space {homogeneous case), but only on time and the magnitude of 
the molecular velocity [isotropic coUisional integral). The same homogeneous 
isotropic case is considered by Truesdell [3] in his famous lectures on natural 
philosophy, as the starting point for investigating the role of time in classical 
thermodynamic systems (which are assumed homogeneous by definition). In 
fact, despite the isotropy of the coUisional integral, the actual time evolutions 
of the distribution function (far from the equilibrium) may be very different. 
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depending on the initial conditions. 



Concerning gas dynamics, focusing on the homogeneous isotropic case, it may 
seem a bit hmiting. For example, an immediate consequence of the isotropic 
symmetry is that all the odd statistical moments are null by definition and 
hence (meaningful) moment equations can be derived for even moments only. 
However it is well known that the decomposition between even moments (pres- 
sure, energy,...) and odd moments (momentum, thermal flux,...) is a key con- 
cept in deriving the fluid dynamic description from the full Boltzmann equa- 
tion, in case of vanishing Knudsen number [B]. In particular, in recovering the 
incompressible limit of the Navier-Stokes equations, the Mach number is as- 
sumed as small as the Knudsen number (diffusive scaling, see [6]) and hence 
the kinetic description collapses in a small neighborhood of the statistical core 
defined by the even moments only. This means that the distribution function 
can be expanded around an equilibrium distribution function, which depends 
on the even moments only. Hence describing properly the manifold defined 
by the even moments is the first basic step for describing the dynamics due 
to small deviations from the local equilibrium. This is the key idea behind 
the derivation of the so-called Lattice Boltzmann Method (LBM) [7]. This is 
the reason why the even moments are sometimes called backbone moments 
of the LBM description [8]. A similar idea holds for the so-called quadrature 
method of moments (QMOM), which is a generic solution method for pop- 
ulation balance models [9]. The common feature between LBM and QMOM 
is that both solve moment systems of equations, which are based on a con- 
traction of the statistical description given by the Boltzmann equation. The 
systematic derivation of moment equations from the Boltzmann equation is 
beyond the purposes of the present work: a detailed review can be found in 
Ref. HO]. 

The interest with regards to the homogeneous isotropic Boltzmann equation 
goes beyond simple dilute gases. In the so-called econophysics [H], a Boltz- 
mann type model is sometimes introduced for studying the distribution of 
wealth in a simple market. The founding idea, dating back to the works of 
Mandelbrot [12], is that the laws of statistical mechanics govern the behavior 
of a huge number of interacting individuals just as well as that of colliding 
particles in a gas container. The classical theory for homogeneous gases is 
easily adapted to the new economic framework: molecules and their velocities 
are replaced by agents and their wealth, and instead of binary collisions, one 
considers trades between two individuals. The goal is to recover the macro- 
scopic distributions of wealth by tuning the microscopic models for the binary 
interaction among the agents. The parameters of the microscopic model can 
be either constant or random quantities. A recent review on this topic can be 
found in Ref. [13] and the references therein. 

Another recent application of the homogeneous isotropic Boltzmann equation 
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is given by opinion formation modeling in quantitative sociology, also called 
sociodynamics or sociophysics [15]. Quantitative sociology has the ambitious 
aim to provide a general strategy, that means a frame of theoretical concepts, 
for designing mathematical models for the quantitative description of a rather 
broad class of collective dynamical phenomena within human society, in par- 
ticular opinion formation. The modeling of opinion dynamics has been treated 
in numerous works, because of its application to politics, to predict the be- 
havior of voters during an election process or the public opinion tendencies 
|15] . Classical kinetic models based on homogeneous isotropic Boltzmann-like 
equations can be derived by prescribing the collision kernel for the micro- 
scopic particle interactions, namely the sociophysical model which prescribes 
the exchange rules for opinion in a binary interaction [13]. 

Since the Boltzmann equation was the starting point for constructing numer- 
ous kinetic equations in many fields of physics, many numerical techniques 
have been proposed to solve it. A complete review of these efforts is clearly 
beyond the purposes of the present work: however a complete discussion can 
be found in the Ref. [16] and the references therein. Despite this wide scenario 
of numerical methods and the constant increase in the computational power, 
solving the Boltzmann equation in practical applications is still challenging 
nowadays. In particular, the most demanding step consists in the evaluation 
of the collisional integral (which is in general five- fold in three dimensions). It 
is possible to distinguish between (a) stochastic and (b) deterministic methods 
in evaluating the collisional integral. In the stochastic methods, like the Monte 
Carlo method, one uses a combination of approximations based on randomly- 
generated variables and a fixed (molecular) velocity grid. On the other hand, 
in deterministic (or direct) methods, one uses only regular lattices in velocity 
space, usually dealing with a larger computational effort in order to achieve 
better accuracy (free of stochastic noise) [TB] . 

The goal of this work is twofold. 

• First of all, this work aims to improve the deterministic numerical method 
proposed by Aristov [16] by (i) reformulating the original problem in terms 
of particle kinetic energy (this allows one to ensure exact particle number 
and energy conservation; momentum is trivially conserved because of the 
isotropic symmetry) and (ii) improving the computation of the relaxation 
rates (making it particularly suitable for dealing with the late dynamics of 
the relaxation towards the equilibrium). 

• Secondly, this work aims to distribute an open-source program (as mini- 
mal as possible) for solving by a deterministic method the homogeneous 
isotropic Boltzmann equation, which can be easily understood and modi- 
fied for dealing with different applications (thermodynamics, econophysics 
and sociodynamics), in order to derive reliable reference solutions (with an 
accuracy which can not be easily obtained by stochastic methods). 
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The paper is organized as follows. First, in Section |2l some theoretical back- 
ground is provided about the homogeneous isotropic Boltzmann equation: in 
particular, the derivation of the homogeneous isotropic case, the energy for- 
mulation, the numerical method, the proposed correction for the relaxation 
rates and the adopted quadrature formulas are discussed. In Sections [3] and 
m an overview of the program structure and a description of the essential 
components are provided. Finally, in Sections [5] and |6l the instructions about 
installation and how to run a test case are provided. 



2 Theoretical background 

2. 1 Boltzmann equation for Maxwell molecules 

Let us consider a dilute gas made of molecules. Let us introduce the probability 
density function, or distribution function, /(t, a;, for the time t G M"*", for the 
position X G M'^, where M'^ is the physical space, and for the molecular velocity 
^ G M|, where ]R| is the velocity space (M^ U ]R| is the phase space). Hence 
the distribution function is defined on the domain {t > 0, a:; G M^, ^ G 
The distribution function allows one to compute the infinitesimal probabil- 
ity to find some molecules in the time interval between t and t + dt, in the 
infinitesimal volume dx & M. around the point x and with a velocity in the 
infinitesimal volume d^ G around the velocity ^, namely f{t,x,^)dtdxd^. 
According to the kinetic theory of gases, the probability density function of 
a dilute gas with elastic binary interactions satisfies the Boltzmann transport 
equation [TP] , namely 

|[ + ^-v./ = g(/,/), (1) 

where the coUisional integral is given by 

Q{fJ)=J J B{g,n)[f{^')f{a-f{^)f{^*)]dnd^., (2) 

M| g n<0 

^* G ]R| is the generic field particle (integration dummy variable) and ci^* 
is its infinitesimal volume in the velocity space; ^'^ G M| are the post- 
collision test and field particle velocities respectively; n G is the unit 
vector along the direction connecting the centers of the two particles during 
the instantaneous collision and versus pointing from particle $, to while 
dn is the infinitesimal solid angle; g' = ^* — ^ is the relative velocity (of the 
field particle with regards to the test particle); finally, B{g, n) is a volumetric 
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particle flux or collision kernel. In the following, we will discuss only the case 
of Maxwell molecules [2] , which highly simplify the expression of the collision 
kernel, namely B{\g ■ n\). Let us assume the following expression 



Bi\g.n\) 



\g- n 



(3) 



where a is the particle radius, c<j is a characteristic mean particle velocity 
(i.e. statistical mean of the particle velocity deviations, which is related to 
the macroscopic sound speed) and 6* is a tunable parameter (natural number, 
i.e. 6* G N). The case 6 = 1 recovers the hard spheres model (popular in 
fluid dynamics), while the case 6 = recovers the constant kernel model 
(which yields constant collision frequency, as commonly done in econophysics 
and sociophysics). In the previous equation, the post-collision test and field 
particle velocities ^' and ^'^ are given by 

^' = ^ + (g.n)n, (4) 
C = ^*-{9- n) n, (5) 

which means that there are many possible outcomes (^',^*) from a given 
pair of incoming (test and field) particle velocities (^,^*), depending on the 
impact direction n obtained by connecting the particle centers during the 
collision. With other words, the generic microscopic collision is defined once 
two additional degree of freedoms are specified (n is a versor). 



2.2 Homogeneous isotropic case 



Let us consider first the homogeneous case (in space). Consequently the proba- 
bility density function becomes f{t, ^) and the homogeneous Boltzmann equa- 
tion becomes 

% = Q{fJ), (6) 



where the coUisional integral is rewritten equivalently as 

Q{fJ)= J J Siq) B{q)[fi^')f{C) - fU)fi^.)]dnd^., (7) 



-oo 



where q = g n, B{q) = B{\q\) for simplicity and S{q) is an auxiliary function 
introduced for simplifying the integration domain (at the price of making more 
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complex the integrand), namely 

S{q) = { (8) 
g > 0. 



Now, let us introduce also the isotropic symmetry of the collision kernel. Be- 
cause of this symmetry, the probability density function is further simplified 
f{t,^), for the time t G M"*" and for the magnitude of the molecular velocity 
= \\$,\\ G M^. In this way, the distribution function allows one to compute the 
infinitesimal probability to find some molecules in the time interval between t 
and t + dt with a velocity magnitude between ^ and ^ + dC,, namely f(t,C,)dtd^. 
Clearly this probability density function can be reformulated in terms of the 
particle kinetic energy E = ^^/2, namely f{t,E). Let us introduce the unit 
vector n^, along the direction and the unit vector Uq along the direction ^, 
namely 



By means of the previous versor, the volume element d^^ can be expressed as 
^Idn^d^^ and consequently 

+00 47r 47r 

QUJ) =111 iff* - ff*) S{q)B{q)edndn.di.. (10) 



It is clear that q is the only parameter potentially dependent on directions n 
and n* and, in general, 

q = $,^-n-i-n = cos (ay) -^cos{a^), (11) 



where is the angle between ^ and n, while Oy is the angle between and n. 
Let us introduce the auxiliary variable x = cos (a^) and y = cos (ay), namely 
q = ^^y — C, X. Let us express the surface elements defined by dn and dn^ in 
Eq. ffTOl) by using ^ as polar axis for dn and n as polar axis for dn^, namely 
dn = sin (a^) da^ df3x and dn^ = sin (ay) day df3y respectively, where /3x and 
f3y are the corresponding azimuthal angles. This yields 

+00 2iT +1 27r +1 

Q{f,f)= J J J J j Wfl- ff*) S{q)B{q)edxd(5,dyd(5ydi.. (12) 

0-10-1 



Taking the square of Eqs. (jH E]) and recalling that q = g ■ n yields 
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i^? =e+q'+2q^x = e{i-x')+e y\ (13) 

{iJf = C + q^-2qi.y = e (1 - y') +ex\ (14) 
Finally, since (/'/^ — //*) does not depend on [i^ and Py, Eq. (fT2l) becomes 
Q{fJ)=N{fJ)-u{f)f, (15) 



where 

+00 +1 +1 

N{fJ)=A7r' J e J J f{e)f{OS{q)B{q)dxdyd^., (16) 

-1 -1 

+ 00 +1 +1 

u{f)=An' J f{^,)e I I S{q)B{q)dxdyd^.. (17) 

-1 -1 



The variables x and y are called collisional parameters (integration dummy 
variables). Let us define Q the domain of integration of the collisional param- 
eters, namely fi=[— 1, 1] x [—1, 1]. It is possible to divide Q in two subregions, 
namely flq>o and fiq<o, defined by g > and g < respectively. Clearly 
f2g>o U f2q<o = ^ and they are separated by the condition q = 0, which is 
the line y = m*x. For any generic point P={xp,yp) G f2g<o, it is possi- 
ble to define another point symmetric with regards to the origin, namely 
P* = {—xp, —yp) G f2g>o. Taking into account that B{q) = B{x,y), it easy to 
prove that B{P) = B{P^) and consequently 

+1+1 +1+1 

j I S{q)B{q)dxdy=- j j B{q)dxdy. (18) 
-1-1 -1-1 



Recalling Eqs. (I13|[T^ . namely 



r = e'(e, e*, y) = ^e{i-x') + ey\ (19) 
e = e:(e, e*, x, y) = y^a^y^y+e^, (20) 

it is easy to prove that ^'(-P) = C,'{P*) and HiP) = ii{P*) and consequently 

+1 +1 +1 +1 

/ / f{e)f{aS{q)B{q)dxdy = - J J f{C)f{OB{q)dxdy. (21) 
-1-1 -1-1 
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Substituting Eq. (EI]) into Eq. and Eq. into Eq. ((IT]) yields 

+00 +1 +1 

N{fJ) = 2n'^a'cl^' J ^ J J fiOfiOl^^y-^^fdxdyd^., (22) 

-1 -1 

+00 +1+1 

uif)= 271' a'cl~' I /(e*)C'/ / \^.y-^xfdxdyd^,. (23) 



-1 -1 



2.3 Energy formulation 



Let us introduce a change of variables in the previous expressions. Let us 
introduce E = ^72, E^ = E' = {i'f /2 and E'^ = {tf /2, namely 

+00 +1+1 

N{f,f) = Ecf I El'' j I f{E')f{E:)\yEl/'-xE'/'fdxdydE,,{24) 



-1 -1 



+1 +1 



u{f) = Ecf J f{E.) El" J J \yEl" - x E'^'f dx dy dE,, (25) 



-1 -1 



where F = 2^^'^^^"7r'a'cs has the dimensions of a volumetric flow rate. Con- 
sequently the collision relations become 

E' = E{l-x') + E,y', (26) 
E: = Ex' + E4l-y'). (27) 

Let us verify the existence of collisional invariants [2] for the previous formu- 
lation. Let us introduce the generic macroscopic quantity $, namely 

+00 

^{t)=ATTV2 J (l){E)fE^/'dE, (28) 



where (f){E) is a generic function of the particle kinetic energy. The macroscopic 
dynamics of the quantity $ can be computed as 

+00 +00 



dt 



J Qif, f) m d$. = 47rv^ j Qif, f) <P{E) E^'' dE, (29) 
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or equivalently 

^ = (0(E), f{E')f{E:) - f{E)m)) , (30) 



where 



, , , +00 +00 +1 +1 



ATTy/2Fc- 



J J J j \yEl'^ ~xE^/^\'^ (i)ip{EE,Y'^dxdydEJE. 



0-1-1 



Clearly the macroscopic dynamics can not depend on the arbitrary labeling 
of the microscopic particles. Hence let us invert E and -E* and, since x = 
cos (ttx) = tiq ■ n and y = cos (ay) = ■ n (see Eqs. Q)), let us invert the 
variables x and y as well. Because of these inversions, the following expression 
holds 

— = (0(E.), fiE')fiE:) - fiE)fiE^)) . (31) 



Next, let us invert the pre- and post-collisional velocities. The coUisional pa- 
rameters expressed by means of the post-collisional velocities become 



It follows immediately that y' ^ E'^ — x' ^TE' = x y/E — y a/E^, which ensures 
that the collisional kernel is unchanged. Equations (l26l 1271 l32l 1331) define the 
transformation {E, E^,,x,y) — {E',E'^,x',y') and they allow one to compute 
the corresponding Jacobian. Its modulus gives the factor by which the trans- 
formation expands or shrinks the infinitesimal volume in the product (0, v^), 
namely 



/ E E 

dx' dy' dE'JE' = J dxdy dE^dE. (34) 



Consequently 

^ = - {<I^{E'), f{E')f{Ei) - f{E)m)) . (35) 
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By using Eqs. (130 1 131 1 135!). it is easy to prove [2] also for the energy formulation 
that 

— = (0(E) + 0(E.) - <P{E') - 0(E:), f{E')f{Ei) - f{E)f{E.)) . (36) 



This means that if the quantity (f){E) is unchanged by the microscopic collision 
(collisional invariant), then the corresponding macroscopic quantity $ is con- 
stant in time (conserved quantity). In particular, let us consider the following 
moments 

%{t) = 47rV2 J f EP+^/^ dE, (37) 



which are obtained by taking (pi^E) = (pp = E^ into Eq. (!28|) . Clearly 0o = 
1 and (pi = E are both invariant during the generic microscopic collision: 
hence, the corresponding macroscopic quantities $o and $i are conserved, 
namely d^o/dt = and d^i/dt = 0. These macroscopic quantities are usually 
formulated in terms of number density 



n 



% = 471^/2 J fE^/'^dE, (3J 



and specific internal energy 



$0 n J ' S^^fE^/^dE- 



The collisional invariants cpQ = 1 and (pi = E (and consequently the conserved 
quantities n and e) are also involved in the definition of the local equilibrium, 
i.e. the distribution function /e such that Qifs, fs) = 0. Let us assume /e = 
exp[— (co(^o + ci0i)], where cq and ci are some proper constants. The collisional 
operator Q{f, f) oc /'/* — //* is consequently null, namely 

QifE, fE) oc exp [-ci{E' + E:)] - exp [-ci{E + E,)] = 0. (40) 



The constants cq and ci can be found by ensuring that Eqs. (138| [39|l are 
satisfied, namely 

^^=(2;^£)^^"K-^)' ^^'^ 
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where Eb = 2e/3. Recalling that the pressure P is defined as one third of the 
stress tensor trace [2], it follows that P = 2/3 ne = uEb- Moreover, recalling 
the ideal gas law, i.e. P = nkBT, where fc^ is the Boltzmann constant and 
T is the temperature, it follows that 2e/3 = Eb = kBT. Introducing the 
specific heat capacity (per mole) at constant volume C^, = e/T, it follows that 
Cy = 3/2 kBi which is correct for monatomic gases considered here. 



2.4 Hierarchy of moment equations 



Sometimes it is more convenient to compute Eq. (1301) in a slightly different 
way, namely 

^ = (0(E), f{E')f{E:)) - {<f>{E), f{E)m)) . (42) 



It has already been shown (in the previous section) that the product (0, ip) 
can be equivalently formulated in terms of the post-coUisional velocities (0, Lf)' 
(since the coUisional kernel is invariant and the infinitesimal volume can 
be transformed by Eq. ( l34l) ). In particular, once the inverse transformation 
{E', El, x\ y') (£", x, y), namely 



E = E'[l-{x'f]+E:{y')\ (43) 

E, = E'{x'f + K[l-{y'fl (44) 

"^^^V^'ii-i^'m^iw ^^^^ 

^^""'iE'ix'Y + Elil-iy'YV ^^^^ 

is used for evaluating 0(-E') = (f){E{E',El,x',y')), the first term in the right 
hand side of Eq. ( 142|) can be rewritten as 

{<P{E)J{E')f{E:)) = {<P{E{E',E:,x',y'))J{E')f{E:)y. (47) 



Omitting the prime symbol in the previous expression allows one to reformu- 
late Eq. (|42]) as 

^ = {<PiE (1 - x') + E^ y') - <P{E), f{E)f{E^)) . (48) 

The previous equation is usually the starting point of the so-called quadra- 
ture method of moments (QMOM), which is a generic solution method for 
population balance models [9]. 
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2. 5 Numerical integration of energy formulation 



Let us assume a maximum value for the test particle kinetic energy E, namely 
Em- Let us divide the interval [0, Em] in M equal parts, with length AE = 
Em/M. Each cell is identified by index 1 < i < M, such that Ei = [i — 
1/2) AE", and the probability distribution function is discretized accordingly, 
namely = f{Ei). As suggested by Ref. [16], this simple discretization (piece- 
wise constant) can be used to compute a numerical approximation Vi of the 
relaxation frequency Vi for the discrete probability distribution function /j, 
namely 



M 

FAEY^f^EfA,, 



(49) 



where 



+1 +1 



Aij = AE-^/^ J J \yE]'^ -xE]'^^ dxdy. 



-1 -1 



(50) 



and 



2(e+3)/2^2^2 



(51) 



The previous expression admits analytical solution, namely 



Ei^' + Ey' 


2+e 


^1/2 _ ^1/2 


2+e 











Ey^Ey\2 + 3e + 9'- 



(52) 



which for 6 = (constant kernel model) yields Aij{0) = 4, while for 6 = 1 
(hard sphere model, consistent with Ref. [I6]) yields 



A,(i) 



1 I 2 Ey^ + 2/3 EiE~^^\ 



E, < E. 



J' 



2 E]'^ + 2/3 > 



(53) 



Similarly, the piecewise discretization can be used to compute a numerical 
approximation Ni of the relaxation frequency Ni for the discrete probability 
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distribution function /«, namely 



M M M 

j=i k=i 1=1 



(54) 



where 



(55) 



and is the compatibility domain (which may also be null). The domain 
Q^j is defined as the locus of points (x, y) & fl such that the post-colhsional 
energies E'{x,y) and El{x,y), defined as 



E\x,y) = Ei{l-x^) + E,y 
Ei{x,y) = E,{l-y^) + E,x 



(56) 
(57) 



are approximated by (piecewise) constants over a small region around the 
point {Ek,Ei). Let us define by E^^ = [k — \)/S.E the lower rounding limit 
and by Ek+ = kAE the higher rounding limit (similarly for Ei^ and Ei^). 
Consequently the pair (x, y) belongs to fifj if {E', E'^) belongs to [Ek^, Ek+\ x 
[£■;_, or equivalently 



Ek- < Ei (1 
El- <Ej{l 



x') + Ej y' < Ek+, 



y 



E,x^<Ei+. 



(58) 
(59) 



Taking into account that E'^ = E'^{E') = Ei + Ej — E', only a segment of the 
function E'^ = E[{E') can (diagonally) fit into the surface element. Hence, in 
order to define fifj, it is enough to solve Eq. (l58l) . which can be reformulated 
as 



^l+ = {{x,y)en:Ei{l 
n_^ = {{x,y) en:Ei{l 



X 



X 



+ E,y^<Ek+}, 
+ Ejy^>Ek-}, 



(60) 
(61) 
(62) 



The regions and fl-cxn are bounded by two hyperbolas and the region 
is the generic intersection between them. This way of defining Qfj is not 
efficient because it requires two different formulas for defining and Q-oo 
respectively. However the problem can be conveniently reformulated, namely 



n_ = {ix,y) en-. Eiii 



X 



E^y^<Ek.}, 



(63) 
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(64) 



In this way, and are defined by the same formula and similarly for the 
integrals defined over them, which can be computed by a unique expression, 
namely 



kl 



(65) 



where 

C{Ek±) 



yE 



1/2 



xE. 



1/2 



dx dy. 



(66) 



In particular, the shape of the domains VL± on the plane (x, y) depend on the 
relative magnitude of the energies Ei, Ej, E^- and Ek+. As it will be discussed 
in the next subsections, six cases are possible, but only three formulas (Ci, C2 
and C3) are required by conveniently switching the arguments, namely 



C{E, 



k± 



Ci{Ei, 


Ej, 


Ek±), 


Ek± < Ei < Ej, 


C2{Ei, 


Ej, 


Ek±), 


Ei < Ek± < Ej, 


CsiEi, 


Ej, 


Ek±), 


Ei < Ej < Ek±, 


Ci{Ej, 


Ei, 


Ek±), 


Ek± < Ej < Ei, 


C2{Ej, 


Ei, 


Ek±), 


Ej < Ek± < Ei, 




Ei, 


Ek±), 


Ej < Ei < Ek±. 



(67) 



Hence, in the following subsections, only the first three cases are discussed. 



2.5.1 Case 1: Ek±<Ei<Ej 

In this case, the domain Q± is defined by 



ai hi 



where a± = yl — Ek±/Ei and b± = y {Ei — Ek±) / Ej. The domain Q± is 
simply made of two strips between two hyperbolas. Taking into account the 
already discussed symmetry of the problem with regards to the origin of the 
plane (x, y), it is possible to save some computations. In particular, considering 
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only the strip such that q{Ei,Ej) < 0, the function Ci{Ei, Ej, Ek±) can be 
expressed as 



+ 1 c±{x) 

Ci{E,,Ej,Eu±) = 2AE-'/^ J J \y E]'^ - x E]'^\' dydx, (69) 

a± -c±(x) 



where 

c±{x) = .J{EiX^-Ei + Ek±)/Ej. (70) 

The previous expression Ci = Ci{6) can be found analytically for particular 
values of 6* G N (a generic expression in terms of 6 was not found). In particular, 
for ^ = (constant kernel model) 

+1 c±{x) 

cm =2 J J {e}^' X - Ef y) dydx = 

a± -c±{x) 

and for 9 = 1 (hard spheres model) 

^i(i) = -Xe^^ J J (^^'^'^ - ^y^y) ^^^^ = 7^v^fe7^-(^2) 

«±-c±(x) « i 



From the previous expressions, if Ek- (< -Efc+) is minimum, i.e. Ek- = 0, then 
Ci(0)=Ci(l) = 0. 



2.5.2 Case 2: Ei < Ek± < Ej 

In this case, the domain Q± is defined by 



y'^ 



where a± = y Ek±/Ei — 1 and &± = y (-Efc± — Ei) /Ej. The domain fi^ is again 
made of two strips between two hyperbolas, but the function C2{Ei, Ej, Ek±) 
can be computed by means of one strip only {q{Ei,Ej) < 0). The integral 
C2 over fl± depends on the coordinates {xi,yi) of the intersections between 
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the previous hyperbolas and yj = ±1. The abscissas of these intersections are 
Xj = ±e± where 




e± = \l+ ' (74) 



In particular, for Ej > Ek±, which is the present case, e± > 1 and consequently 
the intersections (x/, yi) are out of the domain fl±. Hence, for the preset case, 
we can neglect this problem. Consequently 



+1 d-{x) 

C2{E„Ej,E,±) = 2AE~'/^ J I [E^'x-Ey^yY dydx, (75) 

-1 -c±(x) 



where c±{x) is defined by Eq. (1701) and d{x) = x y Ei/Ej. The previous integral 
C2 = 02(9) admits analytical solutions for particular values of ^ e N. In 
particular, for ^ = (constant kernel model) 



«0)^2^H--^.n^-J-^|, (76) 



and for 9 = 1 (hard spheres model) 
3 Ek± — Ej 

3a¥^^e]^' 



^.(l) = 2-^:%7^. (77) 



2.5.3 Case 3: Ei < Ej < Ek± 

In this case, the domain VL^ is defined by 



hi a\ 



where a± = \j Ek±/ Ei — 1 and &± = y {Ek± — E^) /Ej. The domain il± is made 
of a combination of two hyperbolas and the boundaries of ^2: however it is still 
symmetric with regards to the origin and hence the function Cs{Ei, Ej, E^j-) 
can be expressed by means of the subregion with q{Ei,Ej) < 0. Since Ej < 
Ek±i e± < 1 where e± is given by Eq. ( FM|) and consequently the intersections 
{xi,yi) between the previous hyperbolas and yi = ±1 are inside the domain 
fl±. Consequently 
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-e± d{x) 



Cs{Ei, E„ Ek±) = 2 I I ^^1/2 ^ _ ^1/2 ^y^^ ^ 

-1 -1 

2 (e,^/^ X - Ef yf dy dx + 

-e± -c±(x) 

2 Ai?-^/2 I J^K^/^ X - E;/' y)' dydx, 



(79) 



+e± -1 



where c±{x) is given by Eq. flTOj) . (i(x) = x ^jEi/Ej and e± is given by Eq. fl7^ . 
The previous integral C3 = 6*3(6') admits an analytical solution for particular 
values of 6* G N. In particular, for 6' = (constant kernel model) 



C3(0)=4-2 



/ Ei + Ej — E, 



k± 



Ei 



Ej — Ek+ , 

^ l/2pl/2 



E '^E. 



E]'^ + {E, + E, - E,± 
{E^-{E, + Ej-Ek±)f' 



and for 6 = 1 (hard spheres model) 



C3fll 



E. 



1/2 



1 E, 



3< 



?>E]''Ef 



{Ei + Ej — Ek± 



,3/2 



^0) 



Clearly, the previous expressions are always well defined, because the max- 
imum value of Ek_^. (> Ek_) is exactly Ek+ = Ei + Ej, which corresponds 
to = 0. In particular, if Ek+ = Ei + Ej, then 6*3(0) = Aj(0) and 



2.6 Discrete Velocity Model (DVM) and master equation 



As already pointed out, the compatibility domain fifj is defined as the locus 
of points {x,y) G Vt such that, for some given pre-collisional energies {Ei,Ej), 
the post-coUisional energies E\x,y) and E'^{x,y) (see Eqs.fl56l [57|) ) are in 
the neighborhood of the node {Ek, Ei) (coherently with the adopted piecewise 
approximation). Clearly the compatibility domain may also be null. In partic- 
ular, two cases may be distinguished. If the pre-collisional energies are such 
that Ei + Ej < Em, then all the post-collisional energies fit into the adopted 
discretization mesh for the kinetic energy. On the other hand, if Ei + Ej > Em, 
then some post-collisional energies are still physically possible, but they fall 
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outside the discretization mesh (and they should be excluded for consistency). 
Hence purely geometrical considerations yield the following property, namely 



Ei + < Em, (82) 
Ei + Ej > Em-, (83) 

and consequently 



Ei + E^ < Em, (84) 
Ei + Ej > Em. (85) 

In case Ei + Ej < Em, the fact that the equality is exactly satisfied (by 
the discrete numerical operators) is a consequence of the energy formulation, 
which allows one to ensure perfect conservation of particle number and energy 
on a discrete lattice. On the other hand, if Ei + Ej > Em, the pre-collisional 
energies starting from outside of the discretization mesh are automatically 
excluded (even though they are physically possible) and hence also the post- 
coUisional energies falling outside the discretization mesh should be excluded 
as well for consistency. In this way, all the direct and reverse collisions live 
on the same discretization mesh. The latter strategy is advantageous from 
the computational point of view, but it reveals that the adopted numerical 
description only approximates the dynamics due to the collisions with Ei + 
Ej > Em- In case that very accurate simulations are required, it would be 
better to focus on the sub- region [0, £'m/2]. 

Let us define the following matrix 

M M 

4 = EE4'' (86) 

k=l 1=1 



M M 

UU^g = ^' for 

k=l 1=1 
M M 

UU^S<^> for 

k=l 1=1 



for 



M M 

k=l 1=1 
M M 

EY.BS<A„ for 

k=l 1=1 



and consequently Eq. fH^ becomes 



M 

u,^i), = FAEY^ f, Ef A,,. (87) 
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Introducing Qi = Ni — Uifi and taking into account Eqs. (IMl IHTI) yield 

M / M M \ 

= ^ = F 5: E;/' ^ fkfi - 4 . (88) 

7 = 1 \fc=l 1=1 I 



We would like to investigate the (macroscopic) conservation properties of the 
discrete operator Qj. In order to do this, let us rewrite Eq. f l5^ for the dis- 
crete case, which is simplified by the fact that dE'^^ dE' , dE^ and dE are all 
approximated by Ai?, namely 



dx' dy' 



Ei Ej 



\ E'{x,y)Ei{x,y) 



dx dy 



Ek El 



dx dy. 



^9) 



Clearly the last relation is only asymptotically satisfied by the discrete oper- 
ator: hence, even though the microscopic collisions are conservative (in terms 
of mass and kinetic energy), the corresponding macroscopic moments are not 
exactly conserved. The previous relation suggest to multiply and divide Eq. 
( ]88|) by , which allows one to recover the underlaying Discrete Velocity 
Model (DVM) [T7p1 . namely 

§ = ^ E {Ml - m , (90) 

J^i j,k,l=l 



where 



\ / Ej Ea 

ykl _ V ^ T^kl /qi ^1 



The following properties hold [T7] , namely 

= rj^ ^ Tt (92) 



We would like to mention that this kind of models may be affected by the 
problem of (spurious) conservation laws [H]. In this particular case, numerical 
meshes in the velocity/energy space (i.e. lattices) large enough should fix the 

^ In Eq. (|90p . we have adopted a dimensionless F^j, which is different from the 
convention used in Ref. [T7]. However we note that [F\ [E"^/^] [/] = where 

[•] means the physical dimensions. Another difference with regards to Ref. [T7] is 

1/2 

due to the term at the denominator, because of the homogeneous isotropic 

formulation considered here. 
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problem from the practical point of view. Equation (pOi) is sometimes also 
called master equation and F^j is called the matrix of transition frequencies. 

It is possible to correct the matrix of transition frequencies such that it satisfies 
exactly the symmetry properties (DVM correction). There are 24 = 4! possi- 
ble permutations of the four indexes i, j, k and / in the matrix of transition 
frequencies, but only eight permutations ensure the conservation of kinetic 
energy. If two indexes are equal [i = j or k = I), then only four permuta- 
tions (conserving kinetic energy) are possible. Let us define by {T^j} the set 
of transition frequencies obtained by permutations of the indexes {i,j,k,l) 
conserving kinetic energy. The DVM correction is defined as 

V(^,j,fc,/):rgG{rg}, fg = {lf}, (93) 



where the overline means the arithmetic mean of the considered set (sym- 
metrization) . By means of this DVM correction, the following property holds 
(exactly) 

fg = = t%, (94) 



as required by the DVM models [T7]. Consequently it is possible to correct 
the dimensionless frequencies, namely 

gfcf _ f^kl /qc^^ 



M M _ 
k=l 1=1 



which ensure that both particle number and kinetic energy are perfectly con- 
served also at macroscopic level. In the following, the symbols i>i and Qi are 
still used (for keeping the notation as simple as possible), even thought they are 
computed by B^j and Aij instead of B^j and Ay. It is worth the effort to point 
out that, because of the DVM correction, Aij ^ Aij even for Ei + Ej < Em 
(while, under the same conditions, A^j = Aij). 

Ensuring the numerical conservation of conserved hydrodynamic moments is 
also one of the key ideas behind the derivation of the so-called Lattice Boltz- 
mann Method (LBM) [7] (even though mass and momentum only are con- 
served on the smallest lattices). 
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2. 7 Quadrature formulas for computing the moments 

In order to compute the moments defined by Eq. fl37|) . the piecewise constant 
approximation is used. This is consistent with the receipt used for solving the 
colhsional integral Q{f, f) = N{f, f) — i'f- It is worth the effort to point out 
that the property given by Eq. fl94p (and ensured numerically by means of 
the DVM correction) implies the conservation of particle number and energy, 
only if the piecewise constant approximation is used. Hence, even though more 
elaborate quadrature formulas are possible for computing the moments, they 
would spoil the main advantage of the DVM correction, i.e. ensuring that the 
conservation laws are perfectly satisfied. According to the piecewise constant 
approximation, Eq. (ITTj) can be approximated by 

M 

%^~% = 4nV2AE fi E^^^'\ (97) 

i=l 

This way of computing the moments is straightforward, but it produces some 
problems in defining the local equilibrium. Let us suppose to define the local 
discrete equilibrium as (/e)^ = fsiEi), i.e. the local discrete equilibria coincide 
with the nodal values of the continuous function fE defined by Eq. ( 14T]) (for 
some values of n and e). Applying the previous definition yields $o((/£;)i) 7^ n 
and $i((/£;)i) ne, where n and e are defined by continuous integrals in Eq. 
( I38l) and Eq. ( 1391) respectively. This is clearly an effect of the numerical error 
due to the quadrature formula. 

In order to circumvent this problem, let us define the local equilibrium in the 
following way by recursive tuning. For any discrete distribution function /j, 
let us define h = $o(/j) and he = $i(/i). Let us define the 

{fE)i = exp[-{do + 5iEi)], (98) 

where the constants cq and ci are defined such that 

$0 ((/i?)i) = n, <i.i((/^),) =ng. (99) 

By means of this recursive tuning of the local discrete equilibrium, the parti- 
cle number and energy are both constant during the whole relaxation process. 
Eventually, if the continuous distribution function is known as initial condi- 
tion, the assumptions h = $o(/) and he = $i(/) (by Eq. (1371) ) can be used 
instead. 



23 



2.8 Recovering BGK 



It is well known that the coUisional integral of the Boltzmann equation drives 
any initial distribution function towards the local equilibrium [2]. When the 
distribution function is very close to the local equilibrium, the remaining dy- 
namics becomes very slow (on the kinetic time scale) and it can be described 
by the so-called fluid dynamic time scale (which is suitable for describing 
phenomena in the corresponding fluid dynamic regime). Let us search for sim- 
plified expressions of the collisional integral Qi in such regime. The key idea 
is to use the equilibrium distribution function for computing an approxima- 
tion of the relaxation frequency given by Eq. f E7|) . Since we search for an 
approximation of the real relaxation frequency, let us consider Aij (admitting 
analytical expression) instead of Aij in Eq. f l871) . namely z/j ^ {i^E)i where 

_ M 

{UE)^ = FAEY.A.Ey CfE)r (100) 



Consequently, recalling Eq. flHHl) . it is possible to introduce the following ap- 
proximation 



Qi ~ {QB)i = {l^E)i {fE)i - fi 



(101) 



In general, [ue)! still depends on the particle kinetic energy Ei. For fixing the 
ideas, let us consider the Constant Kernel Model - CKM (^ = in Eq. ([3])), 
where A^j (0) = 4 and 

~ W f) W r) 

MO) = 4F (/^), = ^ = (102) 

j=i ny/Z 7rv2 



i.e. the approximated relaxation frequency i^£;(0) is a constant which depends 
on the local number density. A similar procedure can be followed for the Hard 
Sphere Model - HSM (^ = 1 in Eq. ([3])), which also admits an analytical 
expression for (i/£')j(l) involving the error function: see Ref. [2] for details. For 
the present purposes, i.e. the discussion of the numerical results of the test 
case, let us derive the limit of (z/£;)j(l) for high kinetic energies (by considering 
the case E^ > Ej in Eq. (!53l) ). namely 



From the previous limiting case, it is possible to derive the so-called BGK ap- 
proximation [2], where BGK stands for Bhatnagar-Gross-Krook who proposed 
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this simple collisional model. The key idea is to assume a constant relaxation 
frequency (depending on the local number density), namely 



BGK 



^BGK 



{fE)i - fi 



(104) 



In the following, we will assume for simplicity i'bgk = i'^E)i where {i^e)i 
stands for {iyE)i at Ei = AE/2, even though it should be (more precisely) 
i^BGK = limE^o i^e{E). 



3 Overview of the software structure 



In this section, we provide an overview of the HOMISBOLTZ program which was 
developed using Matlab®. The basic idea is to provide a simple illustration of 
the discussed methodology, which can be easily ported to other environments 
(FORTRAN, C++,...). The HOMISBOLTZ program is free software, which can 
be redistributed and/or modified under the terms of the GNU General Public 
License. The HOMISBOLTZ program has been purposely designed in order to 
be minimal, not only with regards to the reduced number of lines (less than 
1,000), but also with regards to the coding style (as simple as possible, hence 
not optimized in terms of execution time). 

A brief flow chart of the program is the following. 
HOMISBOLTZ 

_B_Constant_Kernel_Model(M) , CKM with 6* = in Eq. (j3]) 
B_Constant_Kernel_Model_Creator(M) 

I C_Constant_Kernel_Model(nEi,nEj ,Ek±/AE)= C±, see Eq. ([67]) 

EnergyPerms (t;) , where v = [i,j,k,l]'^ (DVM correction) 

_ Read_B (B ,i ,j ,k,l)= B^j (DVM correction) 
Lwrite_B(B,i, j ,k,l,5*jO (DVM correction) 
_B_Hard_Sphere_Model(M) , HSM with 6' = 1 in Eq. ([3]) 
B_Hard_Sphere_Model_Creator (M) 

I C_Hard_Sphere_Model (nEi , nE j , Ek± /AE) = C± , see Eq. ([67]) 

EnergyPerms (v) , where v = [i,j,k,l]'^ (DVM correction) 

_ Read_B (B,i,j,'k,l)=B^j (DVM correction) 

Lwrite_B(B,i, j ,k,l,5*jO (DVM correction) 
_A_Hard_Sphere_Model(M) , [it may be omitted by Eq. ( l86l) ] 
_ Equilibrium (/i, Ei, AE) = /ij, see Eqs. ([98] [991) 

— Thermodynamics (/j ,AE) = [n,e,kB T] 

— }ivi_Equilihrimi(F, AE,Ei, if E)i, A) = {uE)i, see Eq. ([TOO]) 

— Phi(fi,AE,p)=^p, seeEq. ([97]) 
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Essentially there are two main parts in the HOMISBOLTZ program: (a) comput- 
ing the data structure storing the dimensionless frequencies for the considered 
model, i.e. B(i, j), and (b) the main solution loop (fully explicit and based on 
the forward Euler integration rule). Both parts are described in the following 
sections. 

4 Description of the individual software components 

4-1 Data structure storing the dimensionless frequencies 

The data structure storing the dimensionless frequencies, i.e. B(i,j), is the 
fundamental data structure of the whole program and it is also the most time- 
consuming to be computed. Essentially B(i,j) is the data structure storing 
the dimensionless frequencies B^j (let us suppose that the DVM correction 
applies) for all the GAIN events {E^^Ei) — )■ {Ei,Ej). The dimensionless fre- 
quencies Aij for all the LOSS events {Ei, Ej) — >■ {E^, Ei) can be computed by 
Eq. fl96|) . The matrix B^^ is a four dimensional (sparse) matrix and it is not 
convenient to compute/store it directly. 

Let us introduce a proper labeling for dealing with the sparse matrix B^^. Let 
us define Aj^ the set formed by all the pairs of natural indices (/c, /) such that 
Ek + El = Ei + Ej, with < Ek < Em and < Ei < Em- Let us define 
with Mij the number of elements of the set Aij and let us identify each pair of 
indices by A, namely if 1 < A < Mij then {k{X), /(A)) G Aij. Consequently, Eq. 
(IHH!) (after the DVM correction) can be reformulated by reducing the number 
of nested summations, namely 

Q, = FAEJ2 E]'^ {^^3 - fifj 4) > (105) 

A=l 

and Bij^^''^^^ is simply Bfj for k = k{X) and / = /(A) and it is stored in the 
data structure B(i,j). Hence all the relevant information stored in B(i,j) 
can be labeled by A. 

A brief overview of the data structure B(i, j) is the following. 

B(i,j) 

B ( i , j ) . howmany= M^- 

_B(i,j).k(A)= k{X) 
_B(i,j).l(A)=/(A) 

_B(i,j) .nEkm(A)= Efc(A)-/AE = A;(A) - 1 

_B(i,j) .nEkp(A)= Ek(x)+/AE = k{X) 

— B(i,j) .value (A) = E^/^^'^^^by Eqs. ([651195]). 
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4-2 Main loop 



The main loop of the program aims to compute Qi = Ni — Vi /j at a given time 
step. RecaUing Eq. fllOSp . the operative formulas immediately follow, namely 

M 

N, = FAE'/^ - 1/2)'/' ^^j, (107) 



M 

V, = FAE'/^ Y^ij - l/2)i/V,4- (108) 



The 


previous formulas are implemented straightforwardly in the main 


loop. 


for 


t = . . . ( time ) 






for i = 1:M 






N(i) = 0; 






nu(i) = 0; 






for j = 1:M 






% GAIN term 






Psi(i J) = 0; 






for m = l:B(i,j). howmany 






k = B(i J ).k(m); 






1 = B(i,j).l(m); 






Bijkm = B( i , j ) . value (m) ; 






Psi(i J) = Psi(i ,j)+f(k)*f(l)*Bijkm 






end 






N(i) =N(i) + ... 






F*DeltaE" (3/2)*(j - 1 /2) " ( 1 /2) * Psi ( i , j ) 






% LOSS term 






nu ( i ) = nu ( i ) + . . . 






F*DeltaE" (3/2)*(j -1/2) " (1/2)* f ( j )*A( i 


J ) ; 




end 






Q(i) = N(i)-nu(i)*f(i); 






end 






% Forward Euler integration rule ( explicit ) 






f = f+Deltat.*Q; 




end 
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5 Installation instructions 

The package of the HOMISBOLTZ program consists of five files, namely 

(1) HOMISBOLTZ .m, which is the (single-file) main program (including all the 
subroutines described in the previous flow chart); 

(2) CKM_Structure_B_DVM_nodes_50 .mat, which is the binary data file con- 
taining the data structure B(i, j) in case of the Constant Kernel Model 
(CKM, 6 = in Eq. (^) with M = 50; 

(3) CKM_Structure_B_DVM_nodes_100 .mat, which is the binary data file con- 
taining the data structure B(i, j) in case of the Constant Kernel Model 
(CKM, ^ = in Eq. ([3])) with M = 100; 

(4) HSM_Structure_B_DVM_nodes_50 .mat, which is the binary data file con- 
taining the data structure B(i,j) in case of the Hard Sphere Model 
(HSM, ^ = 1 in Eq. (E])) with M = 50; 

(5) HSM_Structure_B_DVM_nodes_100 .mat, which is the binary data file con- 
taining the data structure B(i,j) in case of the Hard Sphere Model 
(HSM, 6 = 1 in Eq. with M = 100. 

The previous * . mat files are not strictly required. When executed, the pro- 
gram first searches for the binary data file corresponding to the required com- 
bination of collision kernel, discretization resolution (M) and DVM correction 
(ON/OFF). If this binary data file exists, it will be loaded for saving com- 
putational time. Otherwise the data structure B(i,j) will be computed and 
saved as binary data file for future use. Hence the previous binary data files 
are provided as examples. 



6 Test run description 

In this section, a full test case is described. Let us consider a dilute gas made 
of molecules. The interactions among the molecules can be described by means 
of the collision kernel given by Eq. with the following parameters ^ I 

a = l, = 50, 6 = 1 Hard Sphere Model (HSM). (109) 

No BGK-like approximation is adopted, since we want to investigate the full 
nonlinear Boltzmann equation (in the homogeneous isotropic case). 
From the numerical point of view, we can not investigate the full space M"*" for 
the particle kinetic energy. We need to bound our investigations in the range 
[0, Em] and to be sure that the initial conditions well fit into this sub-portion 
of M"*". Actually, in order to achieve better accuracies, as already pointed out 
(see Section [2l6] for details), the whole dynamic phenomenon should fit into 

^ International System of Units (SI) applies. Clearly molecules of 1 m are not 
realistic but this dimension was adopted for simplicity. It is important to point 
out that the characteristic time scale of the relaxation phenomenon r scales as 
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the sub-portion [0, i?j\//2]. Let us consider the following initial condition f{t = 
0,E) = fj{E), where 
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//(/o, Go, Goo) = /o exp 



(110) 



In the considered test case, the values of these parameters are 



Em 



5000 



/o = 5 X 10 



,-4 



Go = 600, 



35. 



(Ill) 



In order to decide the duration of the numerical simulation, we need to in- 
vestigate the characteristic time scale of the relaxation phenomenon r. This 
characteristic time scales as r ~ (nF)^^ ~ {na^Cs)~^. However it may be 
much smaller than that, when the distribution function approaches the local 
equilibrium (fluid dynamic regime). Hence the duration of the phenomenon 
depends also on how far the initial conditions are from the local equilibrium. 
For the present test case (by trials and errors), the duration of the numerical 
simulation was fixed at Tp = 1 x 10~^. 

Finally, the parameters concerning the numerical integration must be spec- 
ified. The range [0, Em] is divided by M = 100 parts and consequently 
AE = Em/M. The time frame Tp is divided by T = 100 part^ and con- 
sequently AT = Tp/T. Since an explicit integration rule is used to solve the 
kinetic equation (namely the forward Euler rule), an upper threshold on the 
discretization time step is expected, namely 



where is a proper constant and 7 is an exponent depending on the mode 
driving the instability (7 = 1 for the advective mode and 7 = 2 for the diffu- 
sive mode). The previous condition is the celebrated Courant-Friedrichs-Lewy 
(CFL) stability condition. The adopted parameters for the considered test 
case satisfy this condition. In the numerical simulations, few non-conserved 
moments are monitored during the relaxation phenomenon, namely $p with 
p G [2 - 9] (see Eq. ^ for details). 

Figure [1] reports the distribution function dynamics from the initial condition 
given by Eq. flllOp . to the local equilibrium given by Eq. (p8 |[99|) . The approach 
to the local equilibrium is initially quite rapid (kinetic stage) and it becomes 
very slow closer to the equilibrium (fluid dynamic stage). It is not so difficult 
to catch the main trend in the dynamics of the distribution function. However 
the formulation in terms of the distribution function may hide some accuracy 
problems in the relaxation of the high-order moments close to the equilibrium. 



^ Temperature is not used directly in the code and this ensures that there is no 
possibility of confusion in the adopted notation. 



AT < k^AE^, 



(112) 
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X 10 




E, particle kinetic energy [m^/s^] 



Fig. 1. (Color online) Distribution function dynamics from the initial condition 
(blue), namely f{t = 0,E) = fi{E) where fi{E) is given by Eq. (jllOp . to the local 
equilibrium (black) , namely given by Eq. (I98| IMl) . 




0.2 0.4 0.6 0.8 1 

Time[s] 



Fig. 2. (Color online) Macroscopic moments dynamics (in time) described by means 
of the relaxation rates Rp given by Eq. (jll3p for p € [2 — 9] . 
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Fig. 3. (Color online) Normalized macroscopic moments dynamics (in time) de- 
scribed by means of the normalized relaxation rates Rp/Rp{t = 0), where Rp is 
given by Eq. (fTTSll for p e [2 - 9]. 

In order to investigate the last point, let us introduce the relaxation rate Rp 
for the macroscopic moment <I>p, namely 



Rr, 



P 



(113) 



where $^ = ^p{fE)- The time evolution of the macroscopic moments with 
p e [2 — 9] is described in Figure |2] by means of Rp and in Figure |3] by means 
of the normalized relaxation rates Rp/Rp{t = 0). Both quantities approach the 
zero value in the late dynamics. The proposed method (and in particular the 
DVM correction and the recursive tuning of the local equihbrium, see Section 
12.71 for details) allows one to catch very precisely the approach to the local 
equilibrium, even by high-order moments. According to the reported results, 
the hard sphere model produces a slower approach to the equilibrium by the 
higher order moments. This point is investigated next. 

In order to check even more precisely the late dynamics of the high-order 
moments, let us introduce a (time-dependent) effectivj^ relaxation frequency 



^ Clearly the definition given by Eq. (11140 leads to an indeterminate form (0/0) for 
t — 7- oo. From the numerical point of view, this may produce some spurious results, 
particularly when the BGK approximation is used. However, this happens when the 
quantity t'p is no more actually relevant. 
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Time[s] ^ 

Fig. 4. (Color online) Time evolution of the normalized effective relaxation frequen- 
cies Vp/uBGR for p G [2 — 9], where Dp is given by Eq. (jll4p and vbgk is the constant 
frequency prescribed by the BGK model (see Section [2.81 for details). 



Vp for the moment p, namely 

(114) 



p 



^p 



where = ^p{Q). In Figure H] the effective relaxation frequencies for p G 
[2 — 9] are normalized by vbgk, where vbgk = {^e)i, {^e)i stands for {h'E)i 
at El = AE/2 and {h'E)i is given by Eq. fllOOl) . The results reported in Figure 
m show that Up < vbgk during the whole dynamics and actually i)p/ ubgk all 
tend to the same asymptotic value (~ 0.665) for t — )■ oo. In order to explain 
such behavior, let us consider the BGK-like approximation given by Eq. f llOip . 
i.e. Qi ~ {QB)i, and let us introduce it in the definition of Up (the subscript i 
has been removed for simplicity), namely 

*p (Je - f) 

The previous approximation allows one to interpret Up as a weighted average 
of the relaxation frequency {i'E)i (valid in the late dynamics) by means of 
the weight (/^ — /)«. The weight {/e — f)i has no definite sign: the ranges 
of Ei where this weight is positive or negative depend on the initial condition 
(both ranges must exist because {fE)i and / have the same number density 
by definition). In particular, the adopted initial condition given by Eq. (11 101) 
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implies [Je ~ f)i < for high kinetic energies (see Figure [T]). Taking into 
account that the relaxation frequency (^'e)^ of the hard sphere model for high 
kinetic energies tends to increase monotonically as ~ y/E'i (according 

to Eq. (1103P ). this leads to a penalization effect in the computation of the 
effective frequency Up. This penalization is larger for higher order moments 
(i.e. it increases with p, as showed in Figure H]) and this explains why the hard 
sphere model produces a slower approach to the equilibrium by the higher 
order moments. 



7 Conclusions 

In this work, some improvements to the deterministic numerical method pro- 
posed by Aristov [16] for the homogeneous isotropic Boltzmann equation are 
discussed. Firstly, the original problem was reformulated in terms of parti- 
cle kinetic energy and this allows one to ensure exact particle number and 
energy conservation during the microscopic collisions (momentum is trivially 
conserved because of the isotropic symmetry). Secondly, the computation of 
the relaxation rates was improved by the DVM correction, which allows one to 
satisfy exactly the macroscopic conservation laws and it is particularly suitable 
for dealing with the late dynamics of the relaxation towards the equilibrium. 
This work aims also to distribute an open-source program (called HOMISBOLTZ), 
which can be easily understood and modified for dealing with different appli- 
cations (thermodynamics, econophysics and sociodynamics), in order to de- 
rive reliable reference solutions (with an accuracy which can not be easily 
obtained by stochastic methods). The HOMISBOLTZ program was developed 
using Matlab@. The basic idea is to provide a simple illustration of the dis- 
cussed methodology, which can be easily ported to other environments (FOR- 
TRAN, C-|--|-,...). The HOMISBOLTZ program is free software, which can be 
redistributed and/or modified under the terms of the GNU General Public 
License. The HOMISBOLTZ program has been purposely designed in order to 
be minimal, not only with regards to the reduced number of lines (less than 
1,000), but also with regards to the coding style (as simple as possible, hence 
not optimized in terms of execution time). 
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